﻿'Süleymaniye Vakfı, Namaz Vakitleri programı, hesaplama işlemleri
'http://www.suleymaniyetakvimi.com/
'Programla ilgili iletişim için : murat@mirasyedi.net
'Hesaplamalarla ilgili iletişim için : 

Imports System.Net
Imports System.Xml
Imports Newtonsoft.Json
Partial Class Hesap
    Inherits System.Web.UI.Page
    Private lat As Double 'Enlem
    Private lng As Double 'Boylam
    Private alt As Double 'Yükseklik

    Private B2 As Double        'Enlem Derece (Enlemin derece değerini, kuzey yarım küre için pozitif, güney yarım küre için negatif olacak şekilde giriniz)
    Private B3 As Double        'Enlem Dakika (Enlemin dakika değerini, kuzey yarım küre için pozitif, güney yarım küre için negatif olacak şekilde giriniz)
    Private B4 As Double        'Boylam Derece (Boylamın derece değerini, doğu yarım küre için pozitif, batı yarım küre için negatif olacak şekilde giriniz)+
    Private B5 As Double        'Boylam Dakika (Boylamın dakika değerini, doğu yarım küre için pozitif, batı yarım küre için negatif olacak şekilde giriniz)
    Private B6 As Integer       'Yıl
    Private B7 As Integer       'ay
    Private B8 As Integer       'gün
    Private B9 As Integer = 2   'Saat bölgesi( Londra'dan saat farkı, ileri ise pozitif, geri ise negatif)
    Private B10 As Integer = 1  'Yaz saati uygulaması (hayır= 0, evet= 1)
    Private B11 As Integer = 0  'Ufuğa göre bölgenin yüksekliği (Ufuk bölgeden yüksekse sıfır, cevabı metre cinsinden yazınız)

    'Sabitler
    Private B13 As Double ' = "00:02:00"    'Öğle namazı kerahat zamanı
    Private C13 As Double '= "06:00:00"    'Öğle namazı kerahat zamanı
    Private B14 As Double                   'Julien Takvim değeri
    Private B15 As Double                   'Asırsal Julien Takvim değeri
    Private B16 As Double                   'Enlem (Ondalık derece)
    Private C16 As Double                   'Enlem (Ondalık derece)
    Private B17 As Double                   'Boylam (Ondalık derece)
    Private B19 As Double                   'Salat Dönencesinin günlük değeri
    Private B20 As Double                   'Seher Vakti (Fecr-i kâzib)
    Private B21 As Double                   'Imsak
    Private B22 As Double                   'Gündoğumu
    Private B23 As Double                   'Öğle
    Private B24 As Double                   'İkindi
    Private B25 As Double                   'Akşam
    Private B26 As Double                   'Yatsı
    Private B27 As Double                   'Yatsı sonu

    'Astronomik değişkenler
    Private AB20 As Double                   'Seher Vakti (Fecr-i kâzib)
    Private AB21 As Double                   'Imsak
    Private AB22 As Double                   'Gündoğumu
    Private AB24 As Double                   'İkindi
    Private AB25 As Double                   'Akşam
    Private AB26 As Double                   'Yatsı
    Private AB27 As Double                   'Yatsı sonu

    Dim objDil As New clsDil
    Protected Sub Page_Load(ByVal sender As Object, ByVal e As EventArgs) Handles Me.Load
        Try
            getValueFromQueryString()

            convertToDegree()

            B13 = 0.00138888888888889
            C13 = 0.25

            Response.Write("<table border='0' width='200'>")
            Response.Write("<tr><td align='center'>")

            GunlukTakvim()

            AylikTakvim()

            MobilUygulama()

            Response.Write("</td></tr></table>")

            'If Request.QueryString("biraylik") = "E" Then

        Catch ex As Exception
            Response.Write(ex.Message)


        End Try

    End Sub

    Private Sub AylikTakvim()
        Dim trh As Date = New Date(B6, B7, B8)

        Response.Write("<div id='tbl4' style='display: none; width: inherit'>")
        Response.Write("<table border='0' width='200' style='font-size:8pt'>")
        Response.Write("<tr>")
        Response.Write(String.Format("<td nowrap='nowrap'><b>{0}</b></td>", objDil.Lang(8)))
        For j As Integer = 20 To 27
            Response.Write(String.Format("<td nowrap='nowrap'><b>{0}</b></td>", objDil.Lang(j)))
        Next
        Response.Write("</tr>")

        For i As Integer = 1 To 30
            B6 = trh.Year
            B7 = trh.Month
            B8 = trh.Day

            doCALC()

            Response.Write("<tr>")
            Response.Write(String.Format("<td align='center'>{0}</td>", trh.ToString("dd/MM/yyyy")))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B20, True)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B21, False)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B22, False)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B23, True)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B24, True)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B25, True)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B26, True)))
            Response.Write(String.Format("<td align='center'>{0}</td>", DecimalToHour(B27, False)))
            Response.Write("</tr>")


            trh = trh.AddDays(1)
        Next
        Response.Write("</table></div>")

    End Sub

    Private Sub GunlukTakvim()
        Dim bln45 As Boolean

        If Math.Abs(B2) >= 45 Then
            bln45 = True
        Else
            bln45 = False
        End If

        doCALC()

        Dim trh As Date = New Date(B6, B7, B8)
        Dim aygunu As Integer = moonage(trh)

        Response.Write("<div id='tbl1' style='width: inherit'>")
        Response.Write("<table border='0' width='100%'>")

        Response.Write("<tr><td>" + objDil.Lang(20) + "</td><td>" + DecimalToHour(B20, True) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(20) + "(ast)</td><td>" + DecimalToHour(AB20, True) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(21) + "</td><td>" + DecimalToHour(B21, False) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(21) + "(ast)</td><td>" + DecimalToHour(AB21, False) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(22) + "</td><td>" + DecimalToHour(B22, False) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(22) + "(ast)</td><td>" + DecimalToHour(AB22, False) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(23) + "</td><td>" + DecimalToHour(B23, True) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(24) + "</td><td>" + DecimalToHour(B24, True) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(24) + "(ast)</td><td>" + DecimalToHour(AB24, True) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(25) + "</td><td>" + DecimalToHour(B25, True) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(25) + "(ast)</td><td>" + DecimalToHour(AB25, True) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(26) + "</td><td>" + DecimalToHour(B26, True) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(26) + "(ast)</td><td>" + DecimalToHour(AB26, True) + "</td></tr>")

        Response.Write("<tr><td>" + objDil.Lang(27) + "</td><td>" + DecimalToHour(B27, False) + "</td></tr>")
        If bln45 Then Response.Write("<tr><td>" + objDil.Lang(27) + "(ast)</td><td>" + DecimalToHour(AB27, False) + "</td></tr>")

        Response.Write("</table></div>")


        Response.Write("<div id='tbl2' style='display: none; width: inherit'>")
        Response.Write("<table border='0' width='100%'>")
        Response.Write("<tr><td>" + objDil.Lang(2) + "</td><td>" + B2.ToString() + ":" + B3.ToString("00") + "</td></tr>")
        Response.Write("<tr><td>" + objDil.Lang(4) + "</td><td>" + B4.ToString() + ":" + B5.ToString("00") + "</td></tr>")
        Response.Write("<tr><td>" + objDil.Lang(11) + "</td><td>" + B11.ToString() + " m</td></tr>")
        Response.Write("<tr><td>" + objDil.Lang(8) + "</td><td>" + B8.ToString() + "/" + B7.ToString() + "/" + B6.ToString() + "</td></tr>")
        Response.Write("<tr><td>" + objDil.Lang(9) + "</td><td>" + B9.ToString() + "</td></tr>")
        Response.Write("<tr><td>" + objDil.Lang(10) + "</td><td>" + IIf(B10 = 0, "Hayır", "Evet") + "</td></tr>")
        Response.Write("</table></div>")

        Response.Write("<div id='tbl3' style='display: none; width: inherit'>")
        Response.Write("<table border='0' width='100%'>")
        Response.Write("<tr><td colspan='2' align='center'><img src='Image/phase" + aygunu.ToString + ".jpg' title='" + aygunu.ToString + "' /><br>" + aygunu.ToString + "</td></tr>")
        Response.Write("</table></div>")

    End Sub

    Private Sub MobilUygulama()
        Response.Write("<div id='tbl5' style='display: none; width: 600px'>")
        Response.WriteFile(Server.MapPath("Mobil/SuleymaniyeTakvimi.htm"))
        Response.Write("</div>")
    End Sub

    Private Sub getDST()
        'If Request.QueryString("yeniyer") <> "E" Then
        '    Return
        'End If

        If getInfoFromGeonames() Then
            Return
        End If

        If getInfoFromEarthTools() Then
            Return
        End If
    End Sub

    Private Function getInfoFromGeonames() As Boolean
        Try
            Dim URI As String = "http://api.geonames.org/timezoneJSON?lat={0}&lng={1}&username=demo"
            URI = String.Format(URI, lat.ToString().Replace(",", "."), lng.ToString().Replace(",", "."))
            Dim wc As New System.Net.WebClient
            Dim strm As System.IO.Stream = wc.OpenRead(URI)
            Dim sr As New System.IO.StreamReader(strm)
            Dim ret As String = sr.ReadToEnd

            Dim values As Newtonsoft.Json.Linq.JObject
            values = JsonConvert.DeserializeObject(ret)

            Dim gmt As String = values("rawOffset").ToString.Trim    '2
            Dim dst As String = values("dstOffset").ToString.Trim    '3

            If gmt.Length > 0 And dst.Length > 0 Then
                If IsNumeric(gmt) And IsNumeric(dst) Then
                    If gmt = dst Then
                        B10 = 0
                    Else
                        B10 = 1
                    End If
                    B9 = CInt(gmt)
                    Return True
                End If
            End If
        Catch ex As Exception

        End Try
        Return False
    End Function

    Private Function getInfoFromEarthTools() As Boolean
        Dim doc As New XmlDocument
        Dim vlr As String
        Try
            Dim url As String = "http://www.earthtools.org/timezone/{0}/{1}"
            url = String.Format(url, lat.ToString().Replace(",", "."), lng.ToString().Replace(",", "."))

            doc.Load(url)

            vlr = doc.SelectSingleNode("//offset").InnerText
            If vlr.Length > 0 And IsNumeric(vlr) Then
                B9 = CInt(vlr)
            End If

            vlr = doc.SelectSingleNode("//dst").InnerText

            If vlr.ToLower = "true" Then
                B10 = 1
            End If
            If vlr.ToLower = "false" Then
                B10 = 0
            End If
            Return True
        Catch ex As Exception

        End Try
        Return False
    End Function

    Private Sub doCALC()
        B14 = calcJD(B6, B7, B8)
        B15 = calcTimeJulianCent(B14)

        B19 = calcSunDeclination(B15)

        If B2 < 45 Then
            B16 = B2 + B3 / 60
        Else
            If B2 < (90 - B19) Then
                If B2 > 0 Then
                    B16 = ((90 - B19) - (B2 + B3 / 60)) * (45 / (45 - B19))
                Else
                    B16 = ((-90 + B19) - (B2 + B3 / 60)) * (45 / (45 - B19))
                End If
            Else
                B16 = (90 - B2 - B3 / 60)
            End If
        End If

        C16 = B2 + B3 / 60
        B17 = B4 + B5 / 60
        B20 = dawn(B16, B17, B6, B7, B8, B9, B10, 18)
        B21 = dawn(B16, B17, B6, B7, B8, B9, B10, 8)

        If 90 - (B2 + B3 / 60) > B19 Then
            B22 = dawn(B16, B17, B6, B7, B8, B9, B10, 0.833 + Math.Sqrt(B11) * 0.0347)
        Else
            B22 = solarnoon(B16, B17, B6, B7, B8, B9, B10) - C13
        End If

        B23 = solarnoon(B16, B17, B6, B7, B8, B9, B10) + B13

        If 90 - (B2 + B3 / 60) > B19 Then
            B25 = dusk(B16, B17, B6, B7, B8, B9, B10, 0.8333 + Math.Sqrt(B11) * 0.0347)
        Else
            B25 = solarnoon(B16, B17, B6, B7, B8, B9, B10) + C13
        End If

        B24 = (B25 + B23) / 2

        B26 = dusk(B16, B17, B6, B7, B8, B9, B10, 9)
        B27 = dusk(B16, B17, B6, B7, B8, B9, B10, 17)

        IkindiNamazi()

        '************************* Astronomik hesaplamalar

        AB20 = dawn(C16, B17, B6, B7, B8, B9, B10, 18)
        AB21 = dawn(C16, B17, B6, B7, B8, B9, B10, 8)
        AB22 = dawn(C16, B17, B6, B7, B8, B9, B10, 0.833 + Math.Sqrt(B11) * 0.0347)
        AB25 = dusk(C16, B17, B6, B7, B8, B9, B10, 0.8333 + Math.Sqrt(B11) * 0.0347)
        AB24 = (AB25 + B23) / 2
        AB26 = dusk(C16, B17, B6, B7, B8, B9, B10, 9)
        AB27 = dusk(C16, B17, B6, B7, B8, B9, B10, 17)
    End Sub

    Private Function DecimalToHour(value As Double, ilave As Boolean) As String
        If value.ToString = "NaN" Then
            Return "??:??"
        End If
        If value > 1 Then
            value = value - 1
        End If
        If ilave Then
            value += (1 / 1440)
        End If
        Dim Deger As Double = value * 24
        Dim Saat As Integer = Math.Truncate(Deger)
        Dim Kalan As Double = (Deger - Saat) * 60
        Dim Dakika As Integer = Math.Truncate(Kalan)
        Dim Saniye As Integer = (Kalan - Dakika) * 60

        'If Saniye > 0 Then
        '    Dakika += 1
        '    If Dakika > 60 Then
        '        Dakika -= 60
        '        Saat += 1
        '        If Saat > 24 Then
        '            Saat -= 24
        '        End If
        '    End If
        'End If
        Return Saat.ToString("00") + ":" + Dakika.ToString("00")
    End Function

    Private Sub getValueFromQueryString()
        Try
            lat = Double.Parse(Request.QueryString("lat"), System.Globalization.CultureInfo.InvariantCulture)
            lng = Double.Parse(Request.QueryString("long"), System.Globalization.CultureInfo.InvariantCulture)
            alt = Double.Parse(Request.QueryString("alt"), System.Globalization.CultureInfo.InvariantCulture)
            B9 = Double.Parse(Request.QueryString("utc"), System.Globalization.CultureInfo.InvariantCulture)
            B10 = Double.Parse(Request.QueryString("dls"), System.Globalization.CultureInfo.InvariantCulture)

            Dim arrTarih() As String = Request.QueryString("trh").Split("-")
            B6 = CInt(arrTarih(0))
            B7 = CInt(arrTarih(1))
            B8 = CInt(arrTarih(2))
            B11 = IIf(alt < 0, 0, alt)

            getDST()

        Catch ex As Exception
            Response.Write(ex.Message)
        End Try

        Try
            If Request.QueryString("yeniyer") <> "E" Then
                Dim cikti As String = ""
                cikti += Request.ServerVariables("REMOTE_ADDR") + "^"
                cikti += Now.ToString("yyyy-MM-dd hh:mm:ss") + "^"
                cikti += "LAT:" + lat.ToString() + "^"
                cikti += "LNG:" + lng.ToString() + "^"
                cikti += "ALT:" + alt.ToString() + "^"
                cikti += Request.ServerVariables("HTTP_USER_AGENT") + "^"

                Log.WriteToLog(cikti)
            End If
        Catch ex As Exception

        End Try
    End Sub

    Private Sub convertToDegree()
        Try
            B2 = Decimal.Truncate(lat)
            B3 = Decimal.Truncate((lat - B2) * 60)

            B4 = Decimal.Truncate(lng)
            B5 = Decimal.Truncate((lng - B4) * 60)
        Catch ex As Exception
            Response.Write(ex.Message)
        End Try
    End Sub

    Private Function radToDeg(angleRad As Double) As Double
        Return (180.0 * angleRad / Math.PI)
    End Function

    Private Function degToRad(angleDeg As Double) As Double
        Return (Math.PI * angleDeg / 180.0)
    End Function

    Private Function calcJD(year As Double, month As Double, day As Double) As Double
        Dim A As Double, B As Double, JD As Double

        If (month <= 2) Then
            year = year - 1
            month = month + 12
        End If

        A = Floor(year / 100, 1)
        B = 2 - A + Floor(A / 4, 1)

        JD = Floor(365.25 * (year + 4716), 1) + _
             Floor(30.6001 * (month + 1), 1) + day + B - 1524.5

        If month = 13 Then
            month = 1
            year = year + 1
        End If
        If month = 14 Then
            month = 2
            year = year + 1
        End If

        Return JD
    End Function

    Private Function calcTimeJulianCent(JD As Double) As Double
        Return (JD - 2451545.0#) / 36525.0#
    End Function

    Private Function calcJDFromJulianCent(t As Double) As Double
        Return t * 36525.0# + 2451545.0#
    End Function

    Private Function calcGeomMeanLongSun(t As Double) As Double
        Dim l0 As Double

        l0 = 280.46646 + t * (36000.76983 + 0.0003032 * t)
        Do
            If (l0 <= 360) And (l0 >= 0) Then Exit Do
            If l0 > 360 Then l0 = l0 - 360
            If l0 < 0 Then l0 = l0 + 360
        Loop

        Return l0
    End Function

    Private Function calcGeomMeanAnomalySun(t As Double) As Double
        Return 357.52911 + t * (35999.05029 - 0.0001537 * t)
    End Function

    Private Function calcEccentricityEarthOrbit(t As Double) As Double
        Return 0.016708634 - t * (0.000042037 + 0.0000001267 * t)
    End Function

    Private Function calcSunEqOfCenter(t As Double) As Double
        Dim m As Double, mrad As Double, sinm As Double, sin2m As Double, sin3m As Double
        Dim c As Double

        m = calcGeomMeanAnomalySun(t)

        mrad = degToRad(m)
        sinm = Math.Sin(mrad)
        sin2m = Math.Sin(mrad + mrad)
        sin3m = Math.Sin(mrad + mrad + mrad)

        c = sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) _
            + sin2m * (0.019993 - 0.000101 * t) + sin3m * 0.000289

        Return c
    End Function

    Private Function calcSunTrueLong(t As Double) As Double
        Dim l0 As Double, c As Double, O As Double

        l0 = calcGeomMeanLongSun(t)
        c = calcSunEqOfCenter(t)

        O = l0 + c
        Return O
    End Function

    Private Function calcSunTrueAnomaly(t As Double) As Double
        Dim m As Double, c As Double, v As Double

        m = calcGeomMeanAnomalySun(t)
        c = calcSunEqOfCenter(t)

        v = m + c
        Return v
    End Function

    Private Function calcSunRadVector(t As Double) As Double
        Dim v As Double, e As Double, R As Double

        v = calcSunTrueAnomaly(t)
        e = calcEccentricityEarthOrbit(t)

        R = (1.000001018 * (1 - e * e)) / (1 + e * Math.Cos(degToRad(v)))
        Return R
    End Function

    Private Function calcSunApparentLong(t As Double) As Double
        Dim O As Double, omega As Double, lambda As Double

        O = calcSunTrueLong(t)

        omega = 125.04 - 1934.136 * t
        lambda = O - 0.00569 - 0.00478 * Math.Sin(degToRad(omega))
        Return lambda
    End Function

    Private Function calcMeanObliquityOfEcliptic(t As Double) As Double
        Dim seconds As Double, e0 As Double

        seconds = 21.448 - t * (46.815 + t * (0.00059 - t * (0.001813)))
        e0 = 23.0# + (26.0# + (seconds / 60.0#)) / 60.0#
        Return e0
    End Function

    Private Function calcObliquityCorrection(t As Double) As Double
        Dim e0 As Double, omega As Double, e As Double

        e0 = calcMeanObliquityOfEcliptic(t)

        omega = 125.04 - 1934.136 * t
        e = e0 + 0.00256 * Math.Cos(degToRad(omega))
        Return e
    End Function

    Private Function calcSunRtAscension(t As Double) As Double
        Dim e As Double, lambda As Double, tananum As Double, tanadenom As Double
        Dim alpha As Double

        e = calcObliquityCorrection(t)
        lambda = calcSunApparentLong(t)

        tananum = (Math.Cos(degToRad(e)) * Math.Sin(degToRad(lambda)))
        tanadenom = (Math.Cos(degToRad(lambda)))


        alpha = radToDeg(Math.Atan2(tanadenom, tananum))

        Return alpha
    End Function

    Private Function calcSunDeclination(t As Double) As Double
        Dim e As Double, lambda As Double, sint As Double, theta As Double

        e = calcObliquityCorrection(t)
        lambda = calcSunApparentLong(t)

        sint = Math.Sin(degToRad(e)) * Math.Sin(degToRad(lambda))
        theta = radToDeg(Math.Asin(sint))
        Return theta
    End Function

    Private Function calcEquationOfTime(t As Double) As Double
        Dim epsilon As Double, l0 As Double, e As Double, m As Double
        Dim y As Double, sin2l0 As Double, sinm As Double
        Dim cos2l0 As Double, sin4l0 As Double, sin2m As Double, Etime As Double

        epsilon = calcObliquityCorrection(t)
        l0 = calcGeomMeanLongSun(t)
        e = calcEccentricityEarthOrbit(t)
        m = calcGeomMeanAnomalySun(t)

        y = Math.Tan(degToRad(epsilon) / 2.0#)
        y = y ^ 2

        sin2l0 = Math.Sin(2.0# * degToRad(l0))
        sinm = Math.Sin(degToRad(m))
        cos2l0 = Math.Cos(2.0# * degToRad(l0))
        sin4l0 = Math.Sin(4.0# * degToRad(l0))
        sin2m = Math.Sin(2.0# * degToRad(m))

        Etime = y * sin2l0 - 2.0# * e * sinm + 4.0# * e * y * sinm * cos2l0 _
                - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m

        Return radToDeg(Etime) * 4.0#
    End Function

    Private Function calcHourAngleDawn(lat As Double, solarDec As Double, solardepression As Double) As Double
        Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

        latRad = degToRad(lat)
        sdRad = degToRad(solarDec)

        HAarg = (Math.Cos(degToRad(90 + solardepression)) / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad))

        HA = (Math.Acos(Math.Cos(degToRad(90 + solardepression)) _
              / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad)))

        Return HA
    End Function

    Private Function calcHourAngleSunrise(lat As Double, solarDec As Double) As Double
        Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

        latRad = degToRad(lat)
        sdRad = degToRad(solarDec)

        HAarg = (Math.Cos(degToRad(90.833)) / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad))

        HA = (Math.Acos(Math.Cos(degToRad(90.833)) _
              / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad)))

        Return HA
    End Function

    Private Function calcHourAngleSunset(lat As Double, solarDec As Double) As Double
        Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

        latRad = degToRad(lat)
        sdRad = degToRad(solarDec)

        HAarg = (Math.Cos(degToRad(90.833)) / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad))

        HA = (Math.Acos(Math.Cos(degToRad(90.833)) _
               / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad)))

        Return -HA
    End Function

    Private Function calcHourAngleDusk(lat As Double, solarDec As Double, solardepression As Double) As Double
        Dim latRad As Double, sdRad As Double, HAarg As Double, HA As Double

        latRad = degToRad(lat)
        sdRad = degToRad(solarDec)

        HAarg = (Math.Cos(degToRad(90 + solardepression)) / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad))

        HA = (Math.Acos(Math.Cos(degToRad(90 + solardepression)) _
               / (Math.Cos(latRad) * Math.Cos(sdRad)) - Math.Tan(latRad) * Math.Tan(sdRad)))

        Return -HA
    End Function

    Private Function calcDawnUTC(JD As Double, latitude As Double, longitude As Double, solardepression As Double) As Double
        Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
        Dim delta As Double, timeDiff As Double, timeUTC As Double
        Dim newt As Double

        t = calcTimeJulianCent(JD)


        eqtime = calcEquationOfTime(t)
        solarDec = calcSunDeclination(t)
        hourangle = calcHourAngleSunrise(latitude, solarDec)

        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta

        timeUTC = 720 + timeDiff - eqtime


        newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0#)
        eqtime = calcEquationOfTime(newt)
        solarDec = calcSunDeclination(newt)
        hourangle = calcHourAngleDawn(latitude, solarDec, solardepression)
        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta
        timeUTC = 720 + timeDiff - eqtime

        Return timeUTC
    End Function

    Private Function calcSunriseUTC(JD As Double, latitude As Double, longitude As Double) As Double
        Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
        Dim delta As Double, timeDiff As Double, timeUTC As Double
        Dim newt As Double

        t = calcTimeJulianCent(JD)


        eqtime = calcEquationOfTime(t)
        solarDec = calcSunDeclination(t)
        hourangle = calcHourAngleSunrise(latitude, solarDec)

        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta

        timeUTC = 720 + timeDiff - eqtime

        newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0#)
        eqtime = calcEquationOfTime(newt)
        solarDec = calcSunDeclination(newt)
        hourangle = calcHourAngleSunrise(latitude, solarDec)
        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta
        timeUTC = 720 + timeDiff - eqtime

        Return timeUTC
    End Function

    Private Function calcSolNoonUTC(t As Double, longitude As Double) As Double
        Dim newt As Double, eqtime As Double, solarNoonDec As Double, solNoonUTC As Double

        newt = calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 + longitude / 360.0#)

        eqtime = calcEquationOfTime(newt)
        solarNoonDec = calcSunDeclination(newt)
        solNoonUTC = 720 + (longitude * 4) - eqtime

        Return solNoonUTC
    End Function

    Private Function calcSunsetUTC(JD As Double, latitude As Double, longitude As Double) As Double
        Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
        Dim delta As Double, timeDiff As Double, timeUTC As Double
        Dim newt As Double

        t = calcTimeJulianCent(JD)


        eqtime = calcEquationOfTime(t)
        solarDec = calcSunDeclination(t)
        hourangle = calcHourAngleSunset(latitude, solarDec)

        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta
        timeUTC = 720 + timeDiff - eqtime



        newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0#)
        eqtime = calcEquationOfTime(newt)
        solarDec = calcSunDeclination(newt)
        hourangle = calcHourAngleSunset(latitude, solarDec)

        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta
        timeUTC = 720 + timeDiff - eqtime


        Return timeUTC

    End Function

    Private Function calcDuskUTC(JD As Double, latitude As Double, longitude As Double, solardepression As Double) As Double
        Dim t As Double, eqtime As Double, solarDec As Double, hourangle As Double
        Dim delta As Double, timeDiff As Double, timeUTC As Double
        Dim newt As Double

        t = calcTimeJulianCent(JD)


        eqtime = calcEquationOfTime(t)
        solarDec = calcSunDeclination(t)
        hourangle = calcHourAngleSunset(latitude, solarDec)

        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta
        timeUTC = 720 + timeDiff - eqtime



        newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC / 1440.0#)
        eqtime = calcEquationOfTime(newt)
        solarDec = calcSunDeclination(newt)
        hourangle = calcHourAngleDusk(latitude, solarDec, solardepression)

        delta = longitude - radToDeg(hourangle)
        timeDiff = 4 * delta
        timeUTC = 720 + timeDiff - eqtime


        Return timeUTC

    End Function

    Private Function dawn(lat As Double, lon As Double, year As Double, month As Double, day As Double, timezone As Double, dlstime As Double, solardepression As Double) As Double
        Dim longitude As Double, latitude As Double, JD As Double
        Dim riseTimeGMT As Double, riseTimeLST As Double


        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8

        JD = calcJD(year, month, day)


        riseTimeGMT = calcDawnUTC(JD, latitude, longitude, solardepression)

        riseTimeLST = riseTimeGMT + (60 * timezone) + (dlstime * 60)

        Return riseTimeLST / 1440
    End Function

    Private Function sunrise(lat As Double, lon As Double, year As Double, month As Double, day As Double, timezone As Double, dlstime As Double) As Double
        Dim longitude As Double, latitude As Double, JD As Double
        Dim riseTimeGMT As Double, riseTimeLST As Double

        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8

        JD = calcJD(year, month, day)


        riseTimeGMT = calcSunriseUTC(JD, latitude, longitude)


        riseTimeLST = riseTimeGMT + (60 * timezone) + (dlstime * 60)

        Return riseTimeLST / 1440
    End Function

    Private Function solarnoon(lat As Double, lon As Double, year As Double, month As Double, day As Double, timezone As Double, dlstime As Double) As Double
        Dim longitude As Double, latitude As Double, JD As Double
        Dim t As Double, newt As Double, eqtime As Double
        Dim solarNoonDec As Double, solNoonUTC As Double


        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8

        JD = calcJD(year, month, day)
        t = calcTimeJulianCent(JD)

        newt = calcTimeJulianCent(calcJDFromJulianCent(t) + 0.5 + longitude / 360.0#)

        eqtime = calcEquationOfTime(newt)
        solarNoonDec = calcSunDeclination(newt)
        solNoonUTC = 720 + (longitude * 4) - eqtime


        solarnoon = solNoonUTC + (60 * timezone) + (dlstime * 60)

        Return solarnoon / 1440
    End Function

    Private Function sunset(lat As Double, lon As Double, year As Double, month As Double, day As Double, timezone As Double, dlstime As Double) As Double
        Dim longitude As Double, latitude As Double, JD As Double
        Dim setTimeGMT As Double, setTimeLST As Double


        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8

        JD = calcJD(year, month, day)


        setTimeGMT = calcSunsetUTC(JD, latitude, longitude)


        setTimeLST = setTimeGMT + (60 * timezone) + (dlstime * 60)

        Return setTimeLST / 1440
    End Function

    Private Function dusk(lat As Double, lon As Double, year As Double, month As Double, day As Double, timezone As Double, dlstime As Double, solardepression As Double) As Double
        Dim longitude As Double, latitude As Double, JD As Double
        Dim setTimeGMT As Double, setTimeLST As Double

        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8

        JD = calcJD(year, month, day)


        setTimeGMT = calcDuskUTC(JD, latitude, longitude, solardepression)


        setTimeLST = setTimeGMT + (60 * timezone) + (dlstime * 60)

        Return setTimeLST / 1440
    End Function

    Private Function solarazimuth(lat As Double, lon As Double, year As Double, month As Double, day As Double, hours As Double, minutes As Double, seconds As Double, timezone As Double, dlstime As Double) As Double

        Dim longitude As Double, latitude As Double
        Dim zone As Double, daySavings As Double
        Dim hh As Double, mm As Double, ss As Double, timenow As Double
        Dim JD As Double, t As Double, R As Double
        Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double
        Dim solarDec As Double, earthRadVec As Double, solarTimeFix As Double
        Dim trueSolarTime As Double, hourangle As Double, harad As Double
        Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double
        Dim azimuth As Double, exoatmElevation As Double
        Dim step1 As Double, step2 As Double, step3 As Double
        Dim refractionCorrection As Double, te As Double, solarzen As Double


        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8

        zone = timezone * -1
        daySavings = dlstime * 60
        hh = hours - (daySavings / 60)
        mm = minutes
        ss = seconds


        timenow = hh + mm / 60 + ss / 3600 + zone

        JD = calcJD(year, month, day)
        t = calcTimeJulianCent(JD + timenow / 24.0#)
        R = calcSunRadVector(t)
        alpha = calcSunRtAscension(t)
        theta = calcSunDeclination(t)
        Etime = calcEquationOfTime(t)

        eqtime = Etime
        solarDec = theta
        earthRadVec = R

        solarTimeFix = eqtime - 4.0# * longitude + 60.0# * zone
        trueSolarTime = hh * 60.0# + mm + ss / 60.0# + solarTimeFix


        Do While (trueSolarTime > 1440)
            trueSolarTime = trueSolarTime - 1440
        Loop

        hourangle = trueSolarTime / 4.0# - 180.0#

        If (hourangle < -180) Then hourangle = hourangle + 360.0#

        harad = degToRad(hourangle)

        csz = Math.Sin(degToRad(latitude)) * _
              Math.Sin(degToRad(solarDec)) + _
              Math.Cos(degToRad(latitude)) * _
              Math.Cos(degToRad(solarDec)) * Math.Cos(harad)

        If (csz > 1.0#) Then
            csz = 1.0#
        ElseIf (csz < -1.0#) Then
            csz = -1.0#
        End If

        zenith = radToDeg(Math.Acos(csz))

        azDenom = (Math.Cos(degToRad(latitude)) * Math.Sin(degToRad(zenith)))

        If (Math.Abs(azDenom) > 0.001) Then
            azRad = ((Math.Sin(degToRad(latitude)) * _
                Math.Cos(degToRad(zenith))) - _
                Math.Sin(degToRad(solarDec))) / azDenom
            If (Math.Abs(azRad) > 1.0#) Then
                If (azRad < 0) Then
                    azRad = -1.0#
                Else
                    azRad = 1.0#
                End If
            End If

            azimuth = 180.0# - radToDeg(Math.Acos(azRad))

            If (hourangle > 0.0#) Then
                azimuth = -azimuth
            End If
        Else
            If (latitude > 0.0#) Then
                azimuth = 180.0#
            Else
                azimuth = 0.0#
            End If
        End If
        If (azimuth < 0.0#) Then
            azimuth = azimuth + 360.0#
        End If

        exoatmElevation = 90.0# - zenith


        '            If (exoatmElevation > 85#) Then
        '                refractionCorrection = 0#
        '            Else
        '                te = Tan(degToRad(exoatmElevation))
        '                If (exoatmElevation > 5#) Then
        '                    refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + _
        '                        0.000086 / (te * te * te * te * te)
        '                ElseIf (exoatmElevation > -0.575) Then
        '                    refractionCorrection = 1735# + exoatmElevation * _
        '                        (-518.2 + exoatmElevation * (103.4 + _
        '                        exoatmElevation * (-12.79 + _
        '                        exoatmElevation * 0.711)))
        '                Else
        '                    refractionCorrection = -20.774 / te
        '                End If
        '                refractionCorrection = refractionCorrection / 3600#
        '            End If


        If (exoatmElevation > 85.0#) Then
            refractionCorrection = 0.0#
        Else
            te = Math.Tan(degToRad(exoatmElevation))
            If (exoatmElevation > 5.0#) Then
                refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + _
                    0.000086 / (te * te * te * te * te)
            ElseIf (exoatmElevation > -0.575) Then
                step1 = (-12.79 + exoatmElevation * 0.711)
                step2 = (103.4 + exoatmElevation * (step1))
                step3 = (-518.2 + exoatmElevation * (step2))
                refractionCorrection = 1735.0# + exoatmElevation * (step3)
            Else
                refractionCorrection = -20.774 / te
            End If
            refractionCorrection = refractionCorrection / 3600.0#
        End If


        solarzen = zenith - refractionCorrection

        '            If (solarZen < 108#) Then
        Return azimuth
        '              solarelevation = 90# - solarZen
        '              If (solarZen < 90#) Then
        '                coszen = Cos(degToRad(solarZen))
        '              Else
        '                coszen = 0#
        '              End If
        '            Else    '// do not report az & el after astro twilight
        '              solarazimuth = -999999
        '              solarelevation = -999999
        '              coszen = -999999
        '            End If

    End Function

    Private Function solarelevation(lat As Double, lon As Double, year As Double, month As Double, day As Double, hours As Double, minutes As Double, seconds As Double, timezone As Double, dlstime As Double) As Double
        Dim longitude As Double, latitude As Double
        Dim zone As Double, daySavings As Double
        Dim hh As Double, mm As Double, ss As Double, timenow As Double
        Dim JD As Double, t As Double, R As Double
        Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double
        Dim solarDec As Double, earthRadVec As Double, solarTimeFix As Double
        Dim trueSolarTime As Double, hourangle As Double, harad As Double
        Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double
        Dim azimuth As Double, exoatmElevation As Double
        Dim step1 As Double, step2 As Double, step3 As Double
        Dim refractionCorrection As Double, te As Double, solarzen As Double


        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8


        zone = timezone * -1
        daySavings = dlstime * 60
        hh = hours - (daySavings / 60)
        mm = minutes
        ss = seconds


        timenow = hh + mm / 60 + ss / 3600 + zone

        JD = calcJD(year, month, day)
        t = calcTimeJulianCent(JD + timenow / 24.0#)
        R = calcSunRadVector(t)
        alpha = calcSunRtAscension(t)
        theta = calcSunDeclination(t)
        Etime = calcEquationOfTime(t)

        eqtime = Etime
        solarDec = theta
        earthRadVec = R

        solarTimeFix = eqtime - 4.0# * longitude + 60.0# * zone
        trueSolarTime = hh * 60.0# + mm + ss / 60.0# + solarTimeFix


        Do While (trueSolarTime > 1440)
            trueSolarTime = trueSolarTime - 1440
        Loop

        hourangle = trueSolarTime / 4.0# - 180.0#

        If (hourangle < -180) Then hourangle = hourangle + 360.0#

        harad = degToRad(hourangle)

        csz = Math.Sin(degToRad(latitude)) * _
              Math.Sin(degToRad(solarDec)) + _
              Math.Cos(degToRad(latitude)) * _
              Math.Cos(degToRad(solarDec)) * Math.Cos(harad)

        If (csz > 1.0#) Then
            csz = 1.0#
        ElseIf (csz < -1.0#) Then
            csz = -1.0#
        End If

        zenith = radToDeg(Math.Acos(csz))

        azDenom = (Math.Cos(degToRad(latitude)) * Math.Sin(degToRad(zenith)))

        If (Math.Abs(azDenom) > 0.001) Then
            azRad = ((Math.Sin(degToRad(latitude)) * _
                Math.Cos(degToRad(zenith))) - _
                Math.Sin(degToRad(solarDec))) / azDenom
            If (Math.Abs(azRad) > 1.0#) Then
                If (azRad < 0) Then
                    azRad = -1.0#
                Else
                    azRad = 1.0#
                End If
            End If

            azimuth = 180.0# - radToDeg(Math.Acos(azRad))

            If (hourangle > 0.0#) Then
                azimuth = -azimuth
            End If
        Else
            If (latitude > 0.0#) Then
                azimuth = 180.0#
            Else
                azimuth = 0.0#
            End If
        End If
        If (azimuth < 0.0#) Then
            azimuth = azimuth + 360.0#
        End If

        exoatmElevation = 90.0# - zenith


        '            If (exoatmElevation > 85#) Then
        '                refractionCorrection = 0#
        '            Else
        '                te = Tan(degToRad(exoatmElevation))
        '                If (exoatmElevation > 5#) Then
        '                    refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + _
        '                        0.000086 / (te * te * te * te * te)
        '                ElseIf (exoatmElevation > -0.575) Then
        '                    refractionCorrection = 1735# + exoatmElevation * _
        '                        (-518.2 + exoatmElevation * (103.4 + _
        '                        exoatmElevation * (-12.79 + _
        '                        exoatmElevation * 0.711)))
        '                Else
        '                    refractionCorrection = -20.774 / te
        '                End If
        '                refractionCorrection = refractionCorrection / 3600#
        '            End If

        If (exoatmElevation > 85.0#) Then
            refractionCorrection = 0.0#
        Else
            te = Math.Tan(degToRad(exoatmElevation))
            If (exoatmElevation > 5.0#) Then
                refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + _
                    0.000086 / (te * te * te * te * te)
            ElseIf (exoatmElevation > -0.575) Then
                step1 = (-12.79 + exoatmElevation * 0.711)
                step2 = (103.4 + exoatmElevation * (step1))
                step3 = (-518.2 + exoatmElevation * (step2))
                refractionCorrection = 1735.0# + exoatmElevation * (step3)
            Else
                refractionCorrection = -20.774 / te
            End If
            refractionCorrection = refractionCorrection / 3600.0#
        End If


        solarzen = zenith - refractionCorrection

        '            If (solarZen < 108#) Then
        '              solarazimuth = azimuth
        Return 90.0# - solarzen
        '              If (solarZen < 90#) Then
        '                coszen = Cos(degToRad(solarZen))
        '              Else
        '                coszen = 0#
        '              End If
        '            Else
        '              solarazimuth = -999999
        '              solarelevation = -999999
        '              coszen = -999999
        '            End If

    End Function

    Private Sub solarposition(lat As Double, lon As Double, year As Double, month As Double, day As Double, hours As Double, minutes As Double, seconds As Double, timezone As Double, dlstime As Double, solarazimuth As Double, solarelevation As Double)
        Dim longitude As Double, latitude As Double
        Dim zone As Double, daySavings As Double
        Dim hh As Double, mm As Double, ss As Double, timenow As Double
        Dim JD As Double, t As Double, R As Double
        Dim alpha As Double, theta As Double, Etime As Double, eqtime As Double
        Dim solarDec As Double, earthRadVec As Double, solarTimeFix As Double
        Dim trueSolarTime As Double, hourangle As Double, harad As Double
        Dim csz As Double, zenith As Double, azDenom As Double, azRad As Double
        Dim azimuth As Double, exoatmElevation As Double
        Dim step1 As Double, step2 As Double, step3 As Double
        Dim refractionCorrection As Double, te As Double, solarzen As Double

        longitude = lon * -1
        latitude = lat
        If (latitude > 89.8) Then latitude = 89.8
        If (latitude < -89.8) Then latitude = -89.8


        zone = timezone * -1
        daySavings = dlstime * 60
        hh = hours - (daySavings / 60)
        mm = minutes
        ss = seconds


        timenow = hh + mm / 60 + ss / 3600 + zone

        JD = calcJD(year, month, day)
        t = calcTimeJulianCent(JD + timenow / 24.0#)
        R = calcSunRadVector(t)
        alpha = calcSunRtAscension(t)
        theta = calcSunDeclination(t)
        Etime = calcEquationOfTime(t)

        eqtime = Etime
        solarDec = theta
        earthRadVec = R

        solarTimeFix = eqtime - 4.0# * longitude + 60.0# * zone
        trueSolarTime = hh * 60.0# + mm + ss / 60.0# + solarTimeFix


        Do While (trueSolarTime > 1440)
            trueSolarTime = trueSolarTime - 1440
        Loop

        hourangle = trueSolarTime / 4.0# - 180.0#

        If (hourangle < -180) Then hourangle = hourangle + 360.0#

        harad = degToRad(hourangle)

        csz = Math.Sin(degToRad(latitude)) * _
              Math.Sin(degToRad(solarDec)) + _
              Math.Cos(degToRad(latitude)) * _
              Math.Cos(degToRad(solarDec)) * Math.Cos(harad)

        If (csz > 1.0#) Then
            csz = 1.0#
        ElseIf (csz < -1.0#) Then
            csz = -1.0#
        End If

        zenith = radToDeg(Math.Acos(csz))

        azDenom = (Math.Cos(degToRad(latitude)) * Math.Sin(degToRad(zenith)))

        If (Math.Abs(azDenom) > 0.001) Then
            azRad = ((Math.Sin(degToRad(latitude)) * _
                Math.Cos(degToRad(zenith))) - _
                Math.Sin(degToRad(solarDec))) / azDenom
            If (Math.Abs(azRad) > 1.0#) Then
                If (azRad < 0) Then
                    azRad = -1.0#
                Else
                    azRad = 1.0#
                End If
            End If

            azimuth = 180.0# - radToDeg(Math.Acos(azRad))

            If (hourangle > 0.0#) Then
                azimuth = -azimuth
            End If
        Else
            If (latitude > 0.0#) Then
                azimuth = 180.0#
            Else
                azimuth = 0.0#
            End If
        End If
        If (azimuth < 0.0#) Then
            azimuth = azimuth + 360.0#
        End If

        exoatmElevation = 90.0# - zenith


        '            If (exoatmElevation > 85#) Then
        '                refractionCorrection = 0#
        '            Else
        '                te = Tan(degToRad(exoatmElevation))
        '                If (exoatmElevation > 5#) Then
        '                    refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + _
        '                        0.000086 / (te * te * te * te * te)
        '                ElseIf (exoatmElevation > -0.575) Then
        '                    refractionCorrection = 1735# + exoatmElevation * _
        '                        (-518.2 + exoatmElevation * (103.4 + _
        '                        exoatmElevation * (-12.79 + _
        '                        exoatmElevation * 0.711)))
        '                Else
        '                    refractionCorrection = -20.774 / te
        '                End If
        '                refractionCorrection = refractionCorrection / 3600#
        '            End If

        If (exoatmElevation > 85.0#) Then
            refractionCorrection = 0.0#
        Else
            te = Math.Tan(degToRad(exoatmElevation))
            If (exoatmElevation > 5.0#) Then
                refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + _
                    0.000086 / (te * te * te * te * te)
            ElseIf (exoatmElevation > -0.575) Then
                step1 = (-12.79 + exoatmElevation * 0.711)
                step2 = (103.4 + exoatmElevation * (step1))
                step3 = (-518.2 + exoatmElevation * (step2))
                refractionCorrection = 1735.0# + exoatmElevation * (step3)
            Else
                refractionCorrection = -20.774 / te
            End If
            refractionCorrection = refractionCorrection / 3600.0#
        End If



        solarzen = zenith - refractionCorrection

        '            If (solarZen < 108#) Then
        solarazimuth = azimuth
        solarelevation = 90.0# - solarzen
        '              If (solarZen < 90#) Then
        '                coszen = Cos(degToRad(solarZen))
        '              Else
        '                coszen = 0#
        '              End If
        '            Else
        '              solarazimuth = -999999
        '              solarelevation = -999999
        '              coszen = -999999
        '            End If

    End Sub

    Public Function Ceiling(ByVal X As Double, Optional ByVal Factor As Double = 1) As Double
        ' X is the value you want to round
        ' is the multiple to which you want to round
        Ceiling = (Int(X / Factor) - (X / Factor - Int(X / Factor) > 0)) * Factor
    End Function

    Public Function Floor(ByVal X As Double, Optional ByVal Factor As Double = 1) As Double
        ' X is the value you want to round
        ' is the multiple to which you want to round
        Floor = Int(X / Factor) * Factor
    End Function

    Function moonage(trh As DateTime) As Double
        Dim j As Double
        Dim ip, ag As Double
        j = juliandate(trh)
        ip = (j + 4.867) / 29.53059
        ip = ip - Math.Floor(ip)
        If (ip < 0.5) Then
            ag = ip * 29.53059 + 29.53059 / 2
        Else
            ag = ip * 29.53059 - 29.53059 / 2
        End If
        ag = Math.Floor(ag) + 1
        Return ag
    End Function

    Function juliandate(trh As DateTime) As Double
        Dim mm, yy, k1, k2, k3, j As Double
        yy = trh.Year - Math.Floor((12 - trh.Month) / 10)
        mm = trh.Month + 9
        If (mm >= 12) Then
            mm = mm - 12
        End If
        k1 = Math.Floor(365.25 * (yy + 4712))
        k2 = Math.Floor(30.600100000000001 * mm + 0.5)
        k3 = Math.Floor(Math.Floor((yy / 100) + 49) * 0.75) - 38
        j = k1 + k2 + trh.Day + 59
        If (j > 2299160) Then
            j = j - k3
        End If
        Return j
    End Function

    Function IkindiNamazi() As Double
        Dim JD As Double = calcJD(B6, B7, B8)
        Dim D As Double = JD - 2451545.0
        Dim a As Double = 357.529 + 0.98560028 * D
        Dim b As Double = 280.459 + 0.98564736 * D
        Dim L As Double = b + 1.915 * Math.Sin(a) + 0.02 * Math.Sin(2 * a)
        Dim JDC As Double = calcTimeJulianCent(JD)
        Dim DC As Double = calcSunDeclination(JDC)

        Dim k As Double = (1 / 15) * Math.Acos((Math.Sin(ArcCot(1 + Math.Tan(L - DC))) - Math.Sin(L) * Math.Sin(DC)) / (Math.Cos(L) * Math.Cos(DC)))
        'K=1/15 arccos⁡((sin⁡(arccot⁡(1+tan⁡(L-DC) ) )-sin⁡(L)×sin⁡(DC))/(cos⁡(L)×cos⁡(DC)))
        'k = 1 / 15 * Math.Acos((Math.Sin(ArcCot(1 + Math.Tan(L - DC))) - Math.Sin(L) * Math.Sin(DC)) / (Math.Cos(L) * Math.Cos(DC)))
        k = Math.Acos((Math.Sin(Math.Atan(1 + Math.Tan(L - DC))) - Math.Sin(L) * Math.Sin(DC)) / (Math.Cos(L) * Math.Cos(DC)))

        Return k
    End Function

    Public Function ArcCot(ByVal x As Double) As Double
        Return Math.PI / 2 - Math.Atan(x)
    End Function

End Class
Public Class Log
    Dim fileLock As New Object()

    Public Shared Sub WriteToLog(text As String)
        Dim fName As String = "Log/" + Now.ToString("yyyyMM") + ".txt"
        fName = System.Web.HttpContext.Current.Server.MapPath(fName)


        SyncLock fName
            Dim tw As System.IO.TextWriter = New System.IO.StreamWriter(fName, True)
            tw.WriteLine(text)
            tw.Close()
        End SyncLock
    End Sub
End Class